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Abstract 

A direct and exact method for calculating the density of states for 
systems with localized potentials is presented. The method is based on 
explicit inversion of the operator E — H. The operator is written in the 
discrete variable representation of the Hamiltonian, and the Toeplitz 
property of the asymptotic part of the obtained infinite matrix is used. 
Thus, the problem is reduced to the inversion of a finite matrix. 
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The evaluation of the density of states has been widely discussed for var- 
ious physical systems in condensed matter and chemical physics 0, EJ. In 
particular, one is interested in the exsistence of resonances, and their posi- 
tion and width. In the field of chemical physics, resonances play a role in 
electron-atom, and atom-molecule scattering processes, autoionization, asso- 
ciative detachment, dissociative atachment, and similar molecular resonant 
reactions [p], |J. Recently, the subject of resonant tunneling in semiconduc- 
tors double-barrier structures has also been the subject of a feverish activity 
[0], due to the technological interest of the properties of such structures, e.g., 
negative differential resistance and bistability in current-voltage response. 

The common methods for finding the position of resonances are the com- 
plex scaling method || and the stabilization method j|, [5]. The complex 
scaling method, pioneered by the works of Aguilar, Balsve and Combes M, 
and Simon J7[ involves the analytic continuation of the energy into the com- 
plex plane using the transformation r — > re te , through diagonalization of 
the scaled (non-Hermitian) Hamiltonian. The following complex eigenval- 
ues correspond to the resonant states, where the real part is the resonance 
energy, and the imaginary part is the life-time. The method was modified 
and extensively used by Certain, Moiseyev and co-workers ||, and has been 
recently applied to three body problems 0. However, this method does not 
give the full density of states as a function of energy. Only an approximated 
form can be obtained through the sum of the Breit-Wigner functions of each 
resonance. 

The stabilization method was developed in quantum chemistry for prob- 
lems involving electron-atom and electron-molecule scattering. The basic 
idea is to repeatedly diagonalize the Hamiltonian in the basis sets of ever 
larger extension L, from what is believed to be the region where the reso- 
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nance wave-function is localized. The result is a stabilization diagram of the 
eigenenergies Ej(L) vs. L (See an example on Figure 1). The characteristic 
plateaus which contain the pattern of avoided crossings between stable (with 
respect to L) and unstable eigenvalues, where the former correspond to the 
eigenvalues representing resonances, and the latter to discretized continuum 
states. Briefly, the explanation of this method is the folowing flnj. The 
wavefunctions for energies near the resonant energy are highly localized, and 
thus are very well described by the finite L 2 states. Therefore, they are sta- 
ble with respect to changes in the range of the basis functions. In contrast, 
eigenvalues far from resonance correspond to extended states, and feel the 
modifications of the basis set. Thus, they will vary as L changes. 

Some methods have been suggested to derive the resonance parameters 
from the stabilization method ITT], |12l ITJ|. More recently, Mandelshtam et 



al (TJJ] have shown how the full density of states can be obtained using some 
kind of averaging over the parameter L. This method has been used for dis- 
sociative photoabsorbsion problems [IS], for the calculation of microcanon- 
ical and canonical rate constants for one- dimensional [1TB] and three-body 



collinear problems [TJJ, and for resonant tunneling-times calculations |T_8|. A 



comparison of the above two methods is given in Ref. |19 | 

However, one would like to have a more direct way to calculate the density 
of states, without applying analytic continuation methods, or (somewhat 
artificial) averaging processes. For instance, the textbook definition of the 
density of states is 

p(E) = Trd(E -H) = — Im TrG(E), (1) 

7T 

where H is the Hamiltonian, and G is the (full) outgoing Green's function 
defined via G(E) = (E + ie — H)~ l . A clear and simple derivation would be 
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the evaluation of the Green's function through the inversion of the matrix 
representing the operator E — H in some L 2 basis set. Unfortunately, this 
naive approach is not directly applicable, since the physical systems in which 
resonant states occur are, in nature, of infinite extent, and accordingly, the 
matrices involved are infinite. A solution for this problem was given by 
Seideman and Miller [pi who applied the method of negative imaginary 



potentials as absorbing boundary conditions |21[ to deal with the infinite 
asymptotes. 

In this Letter, we present an alternative approach for the above prob- 
lem. We manage to invert the infinite matrix E — H without applying to 
any truncations (as done in the stabilization method) or imposing unphysi- 
cal boundary condintions. This is accomplished using the discrete variable 
representation (DVR) |22|, |23], [24]] in which the asyptotic parts are well sepa- 
rated from the interaction region. It has been already recognized |2IJ that the 
asymptotic part of this matrix has a Toeplitz structure This structure 
is here used to reduce the problem into the inversion of a finite matrix whose 
dimesion is proportional to the width of the interaction region. The only 
parameter in this method is the spacing a between successive grid points of 
the DVR. We apply the method to two double barrier problem, with parame- 
ters set corresponding to typical mesoscopic resonant tunneling and chemical 
reaction problems respectively. 

The essence of our method is employing the finite range of the potential 
through the separation between the asymptotic regions and the interaction 
region. In order to keep this separation, a localized basis set is desirable. 
Apparantly, The most appropriate representation from this point of view is 
the DVR, in which the potential operator is diagonal. The representation of 
the kinetic energy part of the Hamiltonian is calculated through the infinite 



4 



order grid point representation of the second derivative, and is (for equally 
spaced grid) of the form 

_ h 2 {-iy-i j tt 2 /3 i = j 

where m is the mass. The corresponding Green's function is thus given 
through 

(G-% = {E- V{x t ))5 l3 - Tij, (3) 

where the points x n = na are the basis grid points. We use in the following 
only the outgoing Green's function, and thus E has to be understood as 
E + ie. 

We now look at the perturbative representation of G 

G = G° + G°VG° + G°VG°VG° + ■■■ = SG°, (4) 

where 

oo 

S=Y / (G°V) n (5) 

n=0 

In the following, we use a block-form in which the vectors (which cor- 
respond to wave-functions) are represented by 3-dimensional super-vectors 
whose first (third) component correspond to the left (right) asymptotic part 
of the vector (and is therefore an infinite-dimensional vector), and its sec- 
ond component correspond to the interaction region (and is thus finite- 
dimensional). Accordingly, the matrices are represented by 3 x 3 super- 
matrices. In this notation, the potential- matrix can be written as (we use 
bold letters for super-matrices elements, to stress that these are matrices 
themselves) 

/ 

V = V 22 I . (6) 

V o o o 
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Similarily, one can write the operator GqV as 

/0 G? 2 V 22 

G°V=\ G° 2 V 22 | . (7) 

V G° 2 V 22 

Since this matrix has a zero (first and third) coloumns, so does all its powers. 
One thus may write 

/I A 

B + I I , (8) 



S 



V C I 

where the I operators are the identity operators of the appropriate order for 
each block. 

It is easy to see from the definition of S, that it satisfies the equation 

S-I = G°VS. (9) 

Explicit multiplication, and comparing term by term, gives the following 
relations 

A = G° 12 V 22 (I + B), (10) 

B = G° 22 V 22 (I + B), (11) 

C = G° 32 V 22 (I + B). (12) 
The solution of these relations is given by 

B = (I-G° 22 V 22 )- 1 G° 22 V 22 , (13) 

and consequently A, C are given by Eqs. fll0l ) ,(|T2| ) . Note that the matrix 
I — G° 22 V 22 is finite and thus its inversion is a simple numerical problem. 
Using (Q), the trace of G is easily obtained. One has 

Tr(G - G°) = Tr(AG 21 ) + Tr(BG 22 ) + Tr(CG° 3 ). (14) 
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The evaluation of the second trace involves simply a finite summation; how- 
ever, the two other traces are infinite sums. We will show now that using the 
explicit form of G° these sums can be reduced to finite ones. 

For this purpose we now calculate the DVR form of G°. By definition, 

G° = (E - T)-\ (15) 

where T stands for the kinetic energy operator. Using DVR, the matrix E—T 
has the structure of a Toeplitz matrix, i.e., {E — T)y = ti_j, where 

tn = \ hH ^Tf , n ■ (16) 
The eigenvalues and eigenvectors are thus given by 

oo *2 2 

\(i)= V t n e inq — E — — (17) 
t#> = -ie***, (18) 

V Z7T 



where g is a continuous index in the region — n < q < n [P7|| . The inverse 
matrix is thus given by 

J-7T AW-' 7tai jo cr + 

where a 2 = —(ak) 2 — ie, and k is the wave number corresponding to the 
energy E. The integrand is highly peaked around q = \a\, and therefore, 
whenever \a\ < it and this peak is inside the integration region, one can 
extend the integration region to infinity, and obtain 

where 9 = ak. The condition 9 < n means that the number of grid points 
per (free) wavelength is > 2. This is a relatively sparse grid with respect to 
those which are used in usual DVR applications. 
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The free particle density of states per unit length is obtained from (f2C|), 
P°(E) = --Ira TiG° = M\ (21) 

71 V 2L TXfl 

Substituting ( p0|) in fll4"|), the infinite sums reduce to geometric serieses, 
and one obtains the formula 

p(E) = -- Y, B mn G^ m +V m (B + I) mn (e^ mM ^^ 

71 m,n=n +l ^ / 1 e 

(22) 

where no (ni) is the first point of the left (right) asymptotic region. Thus, 
after inverting the matrix B + 1, all that one has to do is to sum according 
toEq. (H). 

As an example, we consider here a (symmetric) double barrier structure, 
typical to the problems treated in mesoscopic resonant tunneling problems 
[ |I~8|j . The potential is of the form 

V(x) = Vo(exp(a(a; — x ) 2 ) + exp(a(a; + x ) 2 )) (23) 

The parameters used are Vq = 0.5eV, xo = 50A and a = 4 x 10~ 3 A~ 2 
(corresponding to a quantum well of width ~ 60A, with barriers of width 
~ 40v4). The mass is m* = 0.041m e , corresponding to the Gao.47Ino.53As 
well. Figure 1 shows the full density of states, as calculated form Eq. (|22|) . 
The resulting lineshape is a Lorenzian centered at = 0.101180eK whose 
half-width is T/2 = 1.77127 x 10- 3 eV. Figure 2 presents the same graph 
for the same potential with set of parameters suitable for a chemical physics 
problem, i.e., Vq = 0.5eV, xo = 0.2 A, a = 50 A~ 2 and m = m p . The resulting 
lineshape is a Lorenzian centered at E = 0.3059235eV^ whose half-width is 
T/2 = 3.99496 x 10~ 3 eV. 

In conclusion, a direct and exact method for the calculation of the density 
of states for localized-potential systems was derived. The method evaluates 
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the inverse of the matrix representing the operator E — H in the DVR, 
employing its asymptotic Toeplitz structure. As usual in DVR treatments, 
no integration is needed in constructing the matrix elements. The numerical 
effort needed involves only the inversion of one matrix (for each energy) whose 
dimensionality is proportional to the range of the interaction region. We have 
considered an explicit example of a one-dimensional double barrier structure, 
typical to those considered in the field of resonant tunneling. However, the 
derivation is completely general and can be applied also to more complicated 
(e.g., three-body, three dimensional) problems. 
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Figure Captions: 



Fig. 1: The resonant part of the density of states, i.e., Ap(E) = p(E)—p°(E) 
for the potential fl2"B|). The full curve corresponds to our results, and the 
dotted one is the best Lorentzian fit. 

Fig. 2: Same as Fig. 2 for second set of parameters. 
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Figure 1 - Eisenberg et al 
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Figure 2 - Eisenberg et al 
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